/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Patrick van der Velde
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

using System;

namespace dnAnalytics.LinearAlgebra.Solvers.Preconditioners
{
    /// <summary>
    /// A diagonal preconditioner. The preconditioner uses the inverse
    /// of the matrix diagonal as preconditioning values.
    /// </summary>
    /// <include file='../../../../examples/LinearAlgebra/Solvers/Preconditioners/Diagonal.xml' path='example'/> 
    public sealed class Diagonal : IPreConditioner
    {
        /// <summary>
        /// The inverse of the matrix diagonal.
        /// </summary>
        private double[] m_InverseDiagonals;

        /// <summary>
        /// Returns the decomposed matrix diagonal.
        /// </summary>
        /// <returns>The matrix diagonal.</returns>
        internal SparseMatrix DiagonalEntries()
        {
            SparseMatrix result = new SparseMatrix(m_InverseDiagonals.Length);
            for (int i = 0; i < m_InverseDiagonals.Length; i++)
            {
                result[i, i] = 1 / m_InverseDiagonals[i];
            }
            return result;
        }

        /// <summary>
        /// Initializes the preconditioner and loads the internal data structures.
        /// </summary>
        /// <param name="matrix">
        /// The <see cref="Matrix"/> upon which this preconditioner is based.
        /// </param>
        /// <exception cref="ArgumentNullException">
        /// If <paramref name="matrix"/> is <c>null</c>.
        /// </exception>
        /// <exception cref="MatrixNotSquareException">
        /// If <paramref name="matrix"/> is not a square matrix.
        /// </exception>
        public void Initialize(Matrix matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (matrix.Rows != matrix.Columns)
            {
                throw new MatrixNotSquareException();
            }

            m_InverseDiagonals = new double[matrix.Rows];
            for (int i = 0; i < matrix.Rows; i++)
            {
                m_InverseDiagonals[i] = 1 / matrix[i, i];
            }
        }

        /// <summary>
        /// Approximates the solution to the matrix equation <b>Ax = b</b>.
        /// </summary>
        /// <param name="rhs">The right hand side vector.</param>
        /// <returns>The left hand side vector.</returns>
        public Vector Approximate(Vector rhs)
        {
            if (rhs == null)
            {
                throw new ArgumentNullException("rhs");
            }

            if (m_InverseDiagonals == null)
            {
                throw new MissingMatrixException();
            }

            if (rhs.Count != m_InverseDiagonals.Length)
            {
                throw new NotConformableException();
            }

            Vector result = new DenseVector(rhs.Count);
            Approximate(rhs, result);
            return result;
        }

        /// <summary>
        /// Approximates the solution to the matrix equation <b>Ax = b</b>.
        /// </summary>
        /// <param name="rhs">The right hand side vector.</param>
        /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
        public void Approximate(Vector rhs, Vector lhs)
        {
            if (rhs == null)
            {
                throw new ArgumentNullException("rhs");
            }

            if (lhs == null)
            {
                throw new ArgumentNullException("lhs");
            }

            if (m_InverseDiagonals == null)
            {
                throw new MissingMatrixException();
            }

            if ((lhs.Count != rhs.Count) || (lhs.Count != m_InverseDiagonals.Length))
            {
                throw new NotConformableException();
            }

            for (int i = 0; i < m_InverseDiagonals.Length; i++)
            {
                lhs[i] = rhs[i] * m_InverseDiagonals[i];
            }
        }
    }
}
